Chronic Environmental Perturbation Influences Microbial Community Assembly Patterns

Acute environmental perturbations are reported to induce deterministic microbial community assembly, while it is hypothesized that chronic perturbations promote development of alternative stable states. Such acute or chronic perturbations strongly impact on the pre-adaptation capacity to the perturbation. To determine the importance of the level of microbial pre-adaptation and the community assembly processes following acute or chronic perturbations in the context of hydrocarbon contamination, a model system of pristine and polluted (hydrocarbon-contaminated) sediments was incubated in the absence or presence (discrete or repeated) of hydrocarbon amendment. The community structure of the pristine sediments changed significantly following acute perturbation, with selection of different phylotypes not initially detectable. Conversely, historically polluted sediments maintained the initial community structure, and the historical legacy effect of chronic pollution likely facilitated community stability. An alternative stable state was also reached in the pristine sediments following chronic perturbation, further demonstrating the existence of a legacy effect. Finally, ecosystem functional resilience was demonstrated through occurrence of hydrocarbon degradation by different communities in the tested sites, but the legacy effect of perturbation also strongly influenced the biotic response. This study therefore demonstrates the importance of perturbation chronicity on microbial community assembly processes and reveals ecosystem functional resilience following environmental perturbation.

In this analysis, we used a linear mixed effects model to model community alpha-diversity and include Treatment, Day and an indicator variable HC to denote polluted and pristine sites as fixed effects. The HC indicator variable was constructed by labelling the sites as either pristine or polluted as denoted in figure 2 of the manuscript. We included a three-way interaction between these variables (and all associated two-way interactions) in this model to determine whether alpha diversity changes over time (0 and 28 days), whether this difference is dependent on the treatment (control and phenanthrene) and whether these differences are different between polluted and pristine sites. We included a Site random effect to account for variances associated with site in the model. As data at day 56 was not collected for 'polluted' sites we restricted this initial analysis to Day 0 and Day 28 data. A follow up analysis where we analysed pristine sites across days 0, 28 and 56 is presented in statistic 2.
We then excluded all polluted sites from the original dataframe: env.prist <-env[env$HC == "pristine",] Then we explored several linear models (and checked their assumptions) and the models tested were the following ones: So we have a highly significant two-way interaction between Treatment and Fday with a reported F statistic of 60.094 and an associated P value of P < 0.0001. We concluded that mean Shannon H' does appear to differ over time and this difference is dependent on the treatment.
We then used the emmeans() function from the emmeans package to analyse the results of the model and used the emm object with the pairs() function to produce a Therefore, there is no difference in mean Shannon H' over time for the control group. For the phenanthrene-treated group there is a significant reduction in mean Shannon H' between day 0 and day 56 (difference = -1.2508, P value = 0.0000).

 Statistic 3: Effect of treatment (control or phenanthrene) on evenness (estimated by Pielou's J index) in pristine and polluted sites over time (days 0 and 28):
We applied the same strategy than for statistic 1. We removed one replicate at the North Sea pristine site for Day 0 for consistency between analyses. we explored several linear models (and checked their assumptions) and the models tested were the following ones: env.lm1 <-lm(Pielou ~ Treatment * Fday * HC, data = env.028) env.lme1 <-lme(Pielou ~ Treatment * Fday * HC, random = ~ 1| Sites_HC, data = env.028) env.lme2 <-lme(Pielou ~ Treatment * Fday * HC, random = ~ Fday|Sites_HC,data = env.028) env.lme3 <-lme(Pielou ~ Treatment * Fday * HC, random = ~ Fday| Sites_HC, data = env.028, weights = varIdent(form = ~ 1|Fday)) We selected the model env.lme3 by performing a likelihood ratio test (LRT) and by using AIC. So we have a significant three-way interaction between Treatment, Fday and HC with a reported F statistic of 7.2321 and an associated P value of 0.0076. We concluded that mean Pielou does appear to differ between treatments and this difference is different over time and whether the samples came from a polluted or pristine site.
We then used the emmeans() function from the emmeans package to analyse the results of the model and used the emm object with the pairs() function to produce a Therefore, there is no difference in mean Pielou between control groups at day 0 and day 28 in either pristine or polluted sites (P values = 0.1270 and 0.3892 respectively). In the pristine sites there is a significant reduction in mean Pielou in phenanthrene treated groups between day 0 and day 28 (difference = -2.5020, P value = 0.0130). In comparison, the reduction in Pielou in phenanthrene treated groups between day 0 and day 28 in polluted sites was non-significant (difference = -1.4737, P value = 0.1418).
 Statistic 4: Effect of treatment (control or phenanthrene) on evenness (estimated by Pielou's J index) in pristine sites over time (days 0, 28 and 56): We applied the same strategy than for statistic 2.
We then used the emmeans() function from the emmeans package to analyse the results of the model and used the emm object with the pairs() function to produce a Therefore, there is no difference in mean Pielou over time for the control group. For the phenanthrene-treated group there is a significant reduction in mean Pielou between day 0 and day 56 (difference = -4.6702, P value = 0.0000).

 Statistic 5: Differences between community structure in control and phenanthrene treatments over time
In this analysis we used a PERMANOVA approach to assess community changes between pristine and polluted sites in control and phenanthrene treatments over time. In order to account for multiple measurements at each site we constrained the permutations within sites.
The original OTU data file had five samples removed (rows 57,79,234,298,318) due to low sequence coverage and then rarefied to 9000 for beta diversity measurements using the rrarefy() function from the vegan package. OTUs with a maximum of zero were subsequently removed from these data. A NMDS was then performed on these OTUs using the metaMDSdist() function and then the betadisper() function was used to reduce the original distances to principal components. The scores were then extracted and added to the data in the comm_data.csv data file as a distance variable.
We first excluded day 56 data from the original dataframe (as only available for the pristine sites): env.028 <-env[env$Fday != "56",] env.028$Fday <-factor(env.028$Fday) To keep consistency between this analysis and the alpha-diversity analyses, we removed one replicate at the North Sea pristine site for Day 0.
After importing the OTU table into R and removing the low-read samples, we rarefied all samples to 9000 for beta diversity measurements using the rrarefy() function. Once rarefied, the OTUs with a maximum value of 0 were identified and excluded.

stressplot(all.mds)
Following this, we performed the PERMANOVA on the distance matrix dis using the adonis() function from the vegan package. We also constrained the permutations by site using the strata = argument. We concluded that community similarity is significantly different between control and phenanthrene treated samples and these differences are dependent on day and whether samples were from pristine and polluted sites (P value = 0.001).
We also assessed one of the assumptions of PERMANOVA, which is that the within group variation (dispersion) is similar between groups and we analysed it using the betadisper() function. The dispersions are significantly different between control and phenanthrene treatments. There are also differences between day and whether samples from polluted or pristine sites. While the assumption of similar within group dispersions is not met, it has been reported that adonis() is not highly sensitive to dispersion effects (see ?adonis for more details).

 Statistic 6: Effect of treatment (control or phenanthrene) on community dispersion in pristine and polluted sites over time (days 0 and 28):
We used a similar approach to statistic 1. In this analysis, we used a linear mixed effects model to model community dispersion and include Treatment, Day and an indicator variable HC to denote polluted and pristine sites as fixed effects. The HC indicator variable was constructed by labelling the sites as either pristine or polluted as denoted in figure 2 of the manuscript. We included a three-way interaction between these variables (and all associated two-way interactions) in this model to determine whether dispersion changes over time (0 and 28 days), whether this difference is dependent on the treatment (control and phenanthrene) and whether these differences are different between polluted and pristine sites. We included a Site random effect to account for variances associated with site in the model. As data at day 56 was not collected for 'polluted' sites we restricted this initial analysis to Day 0 and Day 28 data. A follow up analysis where we analysed pristine sites across days 0, 28 and 56 is presented in statistic 7.
To keep consistency between this analysis and the alpha-diversity analyses, we removed one replicate at the North Sea pristine site for Day 0. We also excluded day 56 data from the original dataframe (as only available for the pristine sites) We used the distance matrix using the NMDS transformation of the otunorm object estimated in the previous approach (statistic 5). We subsequently performed the ordination using the betadisper() function and reordered the grouping variables into the two categories (polluted or pristine) dis <-metaMDSdist(otunorm) mod <-betadisper(dis, env$Site, type = c("median","centroid")) Then we explored several linear models (and checked their assumptions) and the models tested were the following ones: env.lm1 <-lm(Distances ~ Treatment * Fday * HC, data = env.028) env.lme1 <-lme(Distances ~ Treatment * Fday * HC, random = ~ 1| Sites_HC, data = env.028) env.lme2 <-lme(Distances ~ Treatment * Fday * HC, random = ~ Fday|Sites_HC, data = env.028) env.lme3 <-lme(Distances ~ Treatment * Fday * HC, random = ~ Fday| Sites_HC, data = env.028, weights = varIdent(form = ~ 1|Fday)) We selected the model env.lme3 by performing a likelihood ratio test (LRT) and by using AIC. So there is a significant three-way interaction between Treatment, Fday and HC with a reported F statistic of 10.9251 and an associated P value of 0.0011. We concluded that mean dispersion does appear to differ between treatments and this difference is different over time and whether the samples came from a polluted or pristine site.
We then used the emmeans() function from the emmeans package to analyse the results of the model and used the emm object with the pairs() function to produce a Therefore, there is no difference in mean dispersion between control groups at day 0 and day 28 in the pristine sites (P value = 0.7037) and the mean dispersion is different in the polluted sites between day 0 and day 28 (P value = 0.0456). In both the pristine and polluted sites, there is no significant difference in mean dispersion in Phenanthrene treated groups between day 0 and day 28 (P value = 0.9568 and 0.1396 respectively). To keep consistency between this analysis and the alpha-diversity analyses, we removed one replicate at the North Sea pristine site for Day 0. We also excluded all polluted sites from the original dataframe.
We used the dispersion index (mod) estimated in the previous approach (statistic 4).
Then we explored several linear models (and checked their assumptions) and the models tested were the following ones: So we have a highly significant two-way interaction between Treatment and Fday with a reported F statistic of 21.5583 and an associated P value of P < 0.0001. We concluded that mean dispersion does appear to differ over time and this difference is dependent on the treatment.
We then used the emmeans() function from the emmeans package to analyse the results of the model and used the emm object with the pairs() function to produce a write.csv(otunorm3, "otunorm3.csv", quote = F) We then selected the 1000 most abundant OTUs in excel and the result is otunorm4.csv We subsequently built a tree on a Biolinux server using the following commands: mafft --auto Lloyd_rarefied.fasta > Lloyd_rarefied_mafft.phy trimal -in Lloyd_rarefied_mafft.phy -out Lloyd_rarefied_trimed.fasta -gt 0.5 iqtree -s Lloyd_mafft_trimed.fasta -alrt 1000 -bb 1000 The tree was then converted in nexus format in figtree and reimported in R phylo <-read.nexus('Lloyd.nex') The phylogenetic tree was then matched to the meanTPH concentration, which was estimated as the ponderated mean of TPH initial (sediment) concentration in function of the abundance of reads within each site.
tph_mean <-read.csv ('OTU_TPHpref.csv') taxa<-data [,1] tph<-data [,2] dim(tph_mean) S17 tph <-tph_mean$taxa The phylogenetic signal was then assessed using the following commands: match.phylo.data(phylo,taxa) phylosignal(tph_mean$tph, phylo, reps = 999, checkdata=TRUE) The mantel correlogram was plotted used the function mantel.correlog from the vegan package: mantel.correlog(D.eco = n, D.geo = x, n.class = 50, cutoff = FALSE, nperm = 999, mult = "bonferroni")  Statistic 9: Phenanthrene degradation over time in pristine versus polluted sediments Due to the destructive sampling approach required for the phenanthrene quantification, phenanthrene degradation was only estimated for polluted sites at day 28 and for pristine sites at day 56. As such, the pristine sites received a double amount of phenanthrene compared to the polluted sites in order to account for the additional time. To take into account the differences in the initial concentration, we calculated the percentage degradation ((start conc -end conc) / start conc) instead of using the final concentration.
Then we explored several linear models (and checked their assumptions) and the models tested were the following ones: end.lme1 <-lme(Percent ~ Polluted, random = ~ 1|Fsite, data = end) end.lme2 <-lme(Percent ~ Polluted, random = ~ 1|Fsite, weights = varIdent(form = ~1|Fsite), data = end) We selected the model end.lme2 by performing a likelihood ratio test (LRT) and by using AIC. We validated the best model (end.lme2) after checking the assumptions: S18 Then, we generated an ANOVA There is a significant difference in mean phenanthrene degradation between pristine and polluted sites (F = 11.76, P value = 0.009).
Finally, we inspected the summary output: